AD-751  566 


ASPECTS  OF  MECHANICAL  BEHAVIOR  OF  ROCK  I 
UNDER  STATIC  AND  CYCLIC  LOADING.  PART  A.  I 
MECHANICAL  BEHAVIOR  OF  ROCK  UNDER  STATIC 
LOADING 

Jesus  Basas,  et  al 
Wisconsin  University 


Prepared  for: 

Bureau  of  Mines 

Advanced  Research  Projects  Agency 


September  1972 


DISTRIBUTED  BY: 


Nation  Technical  Infonutioa  Stnrica 
U.  S.  DEPARTMENT  OF  COMMERCE 

5285  Port  Royal  Road.  Springfield  Va.  22151 


BEST 

AVAILABLE  COPY 


AD 75 1566 


ENGINEERING  EXPERIMENT  STATION 


ASPECTS  OF  MECHANICAL  BEHAVIOR  OF  ROCK 
UNDER  STATIC  AND  CYCLIC  LOADING 

Part  A  Mechanical  Behavior  of  Rock  Under  Static  Loading 

Semi-Annual  Technical  Progress  Report 
September  1972 

Sponsored  by  Advanced  Research  Projects  Agency 
and  monitored  by  U.  S.  Bureau  of  Mines 
under  Contract  No.  HG220041 

ARPA  Order  No.  1579,  Amendment  3 
Program  Code  62701D 


I 


1. 

f  i 

<  « 

•  • 

ASPECTS  OF  MECHANICAL  BEHAVIOR  OF  ROCK  UNDER 
**  STATIC  AND  CYCLIC  LOADING 

-  -  PART  A:  MECHANICAL  BEHAVIOR  OF  ROCK  UNDER  STATIC  LOADING 

Semi-Annual  Technical  Progress  Report 
September  1972 


by 


R.  W.  Heins  (Co- Principal  Investigator 
with  B.  C.  Kaimson) 


Department  of  Metallurgical  &  Mineral  Engineering 

and  the 

Engineering  Exper’ment  Station 
College  of  Engineering 
The  University  of  Wisconsin 
Madison,  Wisconsin  53706 


ARPA  Order  No.  1579,  Amendment  3 
Program  Code  No.  62701D,  Contract  No.  H0220041 
Contract  Period:  March  1972  through  April  1973 
Total  Amount  of  Contract:  $50,  000 
Sponsored  by  ARPA  and  Monitored  by  U.  S.  Bureau  of  Mines 


Disclaimer: 

The  views  and  conclusions  contained  in  this  document  are  those  of 
the  author  and  should  not  be  interpreted  as  necessarily  representing 
the  official  policies,  either  expressed  or  implied,  of  the  Advanced 
Research  Projects  Agency  or  the  U.  S.  Government. 


I 


\ 


uociif.'.t mt  corn koi.  data  .  i:  i  o 
r:,;..*  *-  zzoLzszl  ■- 

_  .  *.  in  i  g«  ■  vi  ivn  I  r  ( i>NKi  li  t  i  iou 

engineering  Experiment  Station  _ Unclassified 

University  of  Wisconsin  /k  ■.*«>.«»  *  “  — — — — 

,  Madison-,  Wisconsin — 53W6- - - - - - 

Aspects  of  Mechanical  Behavior  of  Rock  undei  Static  and  Cyclic  Loading 
_  Part  A  "  Mechanical  Behavior  of  Rock  under  Static  Loading 

4  Ui  ..o i il (Ith.  „r,,,N>.i  JZ iZ (.177,,  J..V.7  - 

— Semi-Annual  Technical  Report 

’•  m  n.mm..  — .... .  .  j  —  .  _ 

Jesus  Basas,  A.  Hayatdavoudi,  Robert  W.  Heins 


a  Hi.fon  t  h  .... 

September  1972 


CWilMCI  O*  CIUM  1  NO - - - - - 

H0220041 

k.  MOJtCI  NO,  .  .  _ 

A  RPA  Order  No.  1579 
«  Amendment  3 

Program  Code  No.  62701D 

rf. 


HnNIlUIIUN  ITAHMlN!  ” 


If*.  TOI  at  tto.  Of  facet  ft#*.  NO.  or  i 


•••  oaiuuAion-i  Kti'OM  mm 

144-C745 


■  s/raw*"  *  N®,k' I4**-  •**'»  »«■—»  «•<.  *•  «•»«#»« 


Distribution  of  this  document  is  unlimited. 


I  II.  IUMI'Umi  «f  Al«v  rout 


I1*-  »Koinvi«iNft  MturiNV  acutrnv 


Advanced  Research  Projects  Agencj 
Washington,  DC  20301 

Brazilian  tests  were  carried  out  to  determine  size-tensile  strength 

Vald"“  Umes,one'  *“»  *•'  C'oud  gray  prtnodiorlte 

laT  thiaM?.  ,  ,hereJS  a  *,fect  0n  *‘r«n*hf  in  Brazilian 

f* f*'  thl®  effect  is  governed  mainly  by  the  position  and  orientation  of  the 
internal  flaws  relative  to  the  loaded  diametral  plane. 

writJr0  Thmen8‘0nal  C°mp"ter  pro*ram  simulating  the  Brazilian  test  was 

based  on  tll  .r  gran  Cmpl°y8  r°ur'sided-  isoparametric  elements  and  is 
based  on  the  same  failure  criteria  used  in  the  first  annual  progress  report 

The  present  program,  however,  is  more  efficient  and  contains8 several  * 

features  not  found  in  the  first  program.  Test  runs  have  proven  that  the 

r/rdiCt  actcurately  the  Progression  of  failure  in  Brazilian  . 
t  and,  to  a  lesser  extent,  the  correlation  of  load  and  displacement  * 

el^pr.°/ram  emPl°yi^  both  two-dimensional  and  three-dimensional 
elements  is  proposed.  The  program  will  be  based  on.  more  realistic 
failure  criteria  and*  will  take  into  account  rock  anisotropy!  ; 


*****  ***«tl>m  In 
X0 


;  •Cm..".  (AM  I  1..*  Kii.:l  l) 

l-  •  v. 


- _ -  - - - - 

i  »•«• 

r 

1  #*.  • 

»  i  ■  *  • 

....  _ 1 

• 

finite  element  analysis 

Brazilian  test 
indirect  tensile  strength 
model  of  rock  failure 
force- displacement  curve 
bounded  distributions  of  strength  and 
stiffness 

finite  element  progressive  failure  model 

1  » 

* 

1 

0 

0 

•  t 

0 

i 

\.  i 

«<•»»  i 

M  < 

Sr  cum  I  y 


preface 


This  report  covers  the 

research  program  entitled,  M  ncipal  Investigator.  The 

Static  entitled. 

PBXaiorof  Rockier  SpraoiecaAo  published  in  a  separate  volume. 
H020041).  Part  B  of  the  projec. 


SUMMARY 


ASPECTS  OF  MECHANICAL  BEHAVIOR  OF  ROCK  UNDER  STATIC  LOADING 

PART  A 


Summary  of  Work  to  Late 

Brazilian  tests  were  carried  out  on  three  rocks  (dacite.  Valders 
limestone,  and  St.  Cloud  gray  granodiorite)  to  determine  size-tensile 
strength  dependence.  Plots  of  tensile  strength  versus  specimen 
dimension  (length  or  diameter)  are  shown  in  Chapter  1.  It  was  concluded 
that  there  is  a  size  effect  on  tensile  strength;  in  Brazilian  test,  this 
effect  is  governed  mainly  by  the  position  and  orientation  of  the  internal 
flaws  relative  to  the  loaded  diametral  plane  rather  than  by  the  extent 
and  number  of  the  flaws. 

A  two-dimensional  computer  program  simulating  the  Brazilian  test 
has  been  completed.  The  program  employs  four-sided,  isoparametric 
elements  and  is  based  on  the  same  failure  criteria  described  in  the  first 
annual  progress  report.  The  present  p  rogram,  however,  is  more 
efficient  and  contains  several  features  not  found  in  the  first  program. 

Test  runs  have  proven  that  the  program  can  predict  accurately  the 
progression  of  failure  in  Brazilian  test  and.  to  a  lesser  extent,  the 

correlation  of  load  and  displacement. 

Development  of  a  program  employing  both  two-dimensional  and 
three-dimensional  elements  has  started.  The  program  will  be  based  on 
more  realistic  failure  criterion  and  will  take  into  account  rock  anisotropy. 
Most  of  the  writing  of  the  program  has  been  done  and  debugging  of  the 
program  is  in  progress. 
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Future  Work 

Development  of  the  combined  two-  or  three-dimensional  program 

will  continue.  Several  examples  will  be  run  to  check  the  correctness 
°*  the  program, 
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CHAPTER  1 


SIZE  EFFECT  ON  BRAZILIAN  TEST 


1  Introduction 

from  lu.  r.  aPP'ying  ralU6S  °f  meChan‘Cal  Pr°Per‘ies  r^S  obtained 
rom  tab  oratory  testa  to  actual  problems,  it  is  essential  for  reasons  of 

safety  and  economy  that  sice  effect,  i,  any,  be  established.  A  design  for 
example,  that  does  not  take  into  account  size  effect  could  be  unsafe  if! 

fact  size  effect  exists.  On  the  other  hand,  a  design  based  on  the  exist¬ 
ence  of  size  effect  could  be  overly  conservative  if  no  such  effect  exists. 

determine  •  T  Cma'de™e  work  has  been  undertaken  to 

.Ton  s‘ze  effect  in  rocks,  the  findings  have  so  far  been  inconclusive 

con  radictory.  In  tests  to  study  size- strength  dependence  it  has 
ee„  „  served  (2)  that,  with  increasing  size,  strength  either  (a)  decreases 
<b>  remains  unchanged,  or  ,c,  increases.  A  very  logical  explanation  of 

'  T  diVergert  effeC‘S  WaS  0ffered  b*  Koifman  ,3).  He  hypo- 
esme  that  size  effect  is  governed  by  two  factors,  namely:  natural 

ernal  taperfections  which  he  called  "volume"  factor  and  "changes  in  the 

rfacelayers.  brought  about  by  mechanical,  physical  or  chemical  action 

y  .  uences  of  the  environment"  which  he  called  "surface"  factor 

"voile"  .C  ard  thi“  W“h  inCrea8ed  SUe-  S‘ren^  ^crease  whe'n  the 
ume  factor  is  dominant,  could  increase  if  the  "surface"  factor  is 

dominant,  or  will  remain  the  same  if  the  two  factors  balance  each  other. 

went  on  to  say  that  under  tensile  stresses  the  "volume"  factor  will  always 

"7“  r  -  ■*-*  --  -» ....... ... : 

and  r®  :?  He  based  hlS  argUment  on  «b*  assumption  that  the  number 
extent  ef  the  internal  Haws  increases  with  size.  Although  this 

assumption  has  a  statistical  basis,  i,  might  not  be  valid  in  actual  situations. 
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Indeed,  in  rock  masses  where  nonhomogeneity  occurs  more  often  that  not 
a  sample  taken  from  a  relatively  defect- free  region  could  easily  contain 
much  fewer  structural  defects  than  smaller  shmples  taken  from  other 
regions  not  quite  as  defect-free.  This  is  particularly  true  for  small-size 
samples  such  as  those  used  in  laboratory  tests.  Furthermore,  in  Brazilian 
tests  the  tensile  strength  could  be  much  more  sensitive  to  the  position 
and  orientation  of  the  structural  defects  relative  to  the  loading  axis  than 
the  number  and  extent  of  these  defects.  Thus,  in  Brazilian  tests  at  least, 
the  possibility  of  strength  increasing  with  size  should  not  be  ruled  out. 

Several  investigators  (2,  3)  have  reported  such  a  size- strength  variation. 

In  the  present  study,  an  attempt  will  be  made  to  correlate  speci¬ 
men  size  and  tensile  strength  as  obtained  from  Brazilian  test.  Three  types  of 
rocks  were  tested.  These  are  Valders  limestone,  St.  Cloud  gray  grano- 
diorite,  and  dacite. 

1. 2  Experimental  Procedure  and  Results 

The  tests  were  carried  out  on  a  MB  Universal  Testing  Machine 
(Fig.  1. 1)  according  to  the  procedure  described  in  the  first  annual  progress 
report  (1).  The  specimens  were  1",  2"  and  3"  in  diameter  and  1/2",  1 '  and 
2"  in  length  (or  thickness),  with  all  nine  possible  combinations  of  length 
and  diameter  represented.  The  specimens  were  tested  in  random  order  and 
without  due  regard  as  to  which  diametral  axis  would  be  loaded.  At  least 
three  samples  of  each  size  were  tested.  A  typical  Brazilian  test  set-up 
is  shown  in  Fig.  1. 2. 

Two  modes  of  testing,  namely,  stress- controlled  mode  and  strain- 
controlled  mode,  were  used.  In  the  stress-controlled  mode,  the  load  is 
applied  at  a  predetermined  constant  rate  of  approximately  100  pounds  per 
second.  In  the  strain- controlled  mode,  the  load  is  applied  at  a  varying 
rate  depending  on  the  lateral  strain  at  the  center  of  the  specimen.  In  the 
latter  mode  the  load  can  be  reduced  faster  than  the  specimen  breaks  thus 
avoiding  the  catastrophic  failure  which  characterizes  stress-controlled  tests. 


I*  injure  1.1 


The  MB  i'imv»-  al  Tt  inn  Machine 


It  is  thus  possible  to  obtain  complete  stress- strain  curves  in  strain- 
controlled  tests. 

To  make  the  two  loading  modes  equivalent  during  the  early  stage 
of  loading,  the  strain-controlled  mode  was  programmed  to  provide  a  maxi¬ 
mum  strain  of  10,  000  micro-inch  per  inch  in  800  seconds.  This  rate,  it 
was  estimated,  is  approximately  equal  to  the  100  lbs.  /  sec.  rate  used  in 
the  stress-controlled  mode. 

Plots  of  tensile  strength  (o^)  versus  specimen  dimensions  for  all 
rocks  are  shown  in  Figs.  1.  3  through  1. 14.  The  curves  were  plotted  on  the 
basis  of  the  equation, 

0  2  0  3 
«t  =  OjL  2  D 


in  which  L  is  length;  D  is  diameter;  and  0^,  O^,  and  0^  are  constant  para¬ 
meters.  The  values  of  the  constant  parameters  corresponding  to  the 
condition  of  "best  fit"  can  be  obtained  by  means  of  a  statistical  procedure 
called  regression  analysis  (5). 


1. 3  Discussion  of  Results 

All  kinds  of  strength- size  variation  are  shown  in  the  plots.  In  the 
harder  rocks  (St.  Cloud  gray  granodiorite  and  Valders  limestone  in  some 
cases)  there  appears  to  be  a  definite  correlation  between  size  and  strength. 
Some  agree  with  Koifman's  prediction.  The  apparent  absence  of  definite 
pattern  of  the  size-strength  relationship  in  the  other  plots  could  be  attributed 
to  one  or  a  combination  of  the  following; 

(a)  The  size  differences  between  the  specimens  were  not  large 
enough. 

(b)  Not  enough  samples  were  tested  for  certain  sizes  parti¬ 
cularly  in  strain-controlled  tests. 

(c)  Valders  limestone  and  dacite  are  not  nearly  homogenous 
and  isotropic  as  first  thought.  The  tensile  strength  waa 
therefore,  affected  more  by  the  position  and  orientation 
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of  the  internal  structural  defects  with  respect  to  the 
loaded  diametral  plane  than  by  the  number  and  extent 
of  said  defects.  This  also  explains  the  scattering  of 
the  data  points. 

1. 4  Conclusions 

There  is  definitely  a  size  effect  on  tensile  strength.  In 
Brazilian  tests,  the  size  effect  will  be  governed  mainly  by  the  position 
and  orientation  of  the  internal  defects  relative  to  the  loaded  diametral 
plane  rather  than  by  the  extent  and  number  of  the  defects.  For  this  reason, 
the  Brazilian  test  is  not  a  good  basis  for  studying  size-tensile  strength 
dependence  unless  care  is  taken  to  consistently  load  the  specimens  elong 
the  same  diametral  plane.  The  splitting  test  described  by  Koifman  (6) 
appears  to  be  a  better  alternative. 

The  size  difference  between  the  specimens  used  in  this  study 
was  not  large  enough  to  predict  a  definite  pattern  of  size-strength  relation¬ 
ship  especially  in  the  case  of  the  soft  rocks.  Future  tests  should  include 
larger  specimens. 
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CHAPTER  2 

FINITE  ELEMENT  SIMULATION  OF  BRAZILIAN  TEST _ 

2.  1  Introduction 

A  theoretical  version  of  the  Brazilian  test  employing  the  finite  element 
technique  is  described  in  this  chapter.  In  connection  with  the  proposed 
method,  development  of  a  highly  efficient  computer  program  which  will  be 
based  on  both  three-dimensional  and  two-dimensional  elements  and  which 
will  account  for  material  anisotropy  as  well  as  nonhomogeneity  is  now  under¬ 
way.  A  two-dimensional  program  which  takes  into  account  nonhomogeneity 
but  not  anisotropy  has  already  been  developed  and  will  be  described  in  the 

next  chapter. 

In  the  treatment  of  nonhomogeneous  problems,  the  basic  idea  suggested 
in  the  first  annual  report  will  be  used;  that  is,  elastic  properties  will  be 
assigned  randomly  to  each  element  by  means  of  a  random  number  generating 
routine.  Other  failure  criteria  not  touched  in  the  first  annual  report  will  be 
investigated.  It  has  been  experimentally  demonstrated  (1)  that  although 
rock  is  essentially  a  brittle  material,  it  attains  an  unusually  high  degree  of 
ductility  when  subjected  to  high  confining  pressures.  This  phenomenon, 
however,  is  vaguely  defined  in  literature  and  will,  therefore,  be  taken  into 
account  only  approximately  when  considering  post-failure  behavior  of  elements. 

In  writing  this  report,  it  is  assumed  that  the  reader  is  familiar  with 
the  basic  principles  of  matrix  algebra  and  the  finite  element  method.  No 
attempt  will  be  made  to  rederive  equations  which  have  already  been  derived 
in  previous  publications. 

2.  2  Formulation  of  Essential  Matrices  for  Three-Dimensional  Elements, 

In  the  finite  element  method  of  analysis,  the  whole  structural  system 
is  idealized  as  an  assemblage  of  elements  which  are  connected  to  one  another 
only  at  a  discrete  number  of  nodal  points.  The  nodal  displacements  are  the 
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basic  unknown  quantities  of  the  method  upon  which  the  displacement  pattern 
and  therefore  the  stresses  within  the  boundaries  of  the  element  depend.  To 
facilitate  the  calculation  of  the  element  stiffness  properties,  a  set  of  dis¬ 
placement  functions,  usually  polynomials,  are  assumed.  These  functions 
uniquely  define  the  deformations  allowed  within  an  element  in  terms  of  the 
nodal  displacements. 

In  this  study,  the  Brazilian  test  specimen  is  divided  by  imaginary 
annular  surfaces  and  radial  planes  into  elements  of  the  kind  shown  in  Fig.  2.  2. 
The  comer  points  of  the  elements  are  designated  as  the  nodal  points.  Eight- 
term  linear  polynomials  are  used  to  represent  the  radial  (u),  circumferential 
(v),  and  axial  (w)  displacements  within  an  element.  Thus 
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(2.  1) 


where  a^,  a2> . .  a24  are  the  constant  coefficients  of  the  polynomials 

and  (r\  0',  z')  are  local  dimensionless  cylindrical  coordinates  which  range  in 
value  from  -1  to  +1  within  an  element  The  global  cylindrical  coordinates 
(r,  0,  z)  are  related  to  the  local  coordinates  as  follows: 


r  =  r  +  r  r' 
o  s 

0  =  0  4-0  0' 
o  s 

z  =  z  +  z  z' 
o  s 


(2.  2) 


wl'*rc  (rQ,  0q,  z q)  are  the  global  coordinates  of  the  origin  of  the  local  coordin¬ 
ate  axes  and  2rg,  20g  and  2zg  are  the  side  dimensions  of  die  element  (see 
Fig.  2.  3). 


FIGURE  2.  2  Element  Models  With  Local  Node  Numbers  Indicated 
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To  obtain  the  displacements  at  any  node  i,  one  merely  substitutes 
into  eq.  (2.  1)  the  appropriate  nodal  coordinates,  that  is, 


u,  *  ®i + vi + a38i + v;8; + 


vi  *  a9  +  aiori  +  an8i  +  ai2ri8i  + 


(2.3) 


wi  =  a17  +  a18ri  +  a199i  +  a20ri8i  + 


where  rj,  0J,  zj  are  equal  to  +1  or  -1.  Thus,  if  (x)  denotes  the  nodal  dis¬ 
placement  vector  (Uj  Vj  w ^  u2  v2  - vQ  Wg)  and  la}  the  constant  co¬ 

efficient  array,  then 

lx}  =  [c]  la} 

The  matrix  [c]  is  shown  in  Table  2.  1.  It  can  easily  be  verified  that 

,-l 


[c] 


1  r  l  T 
¥  c 


From  the  equations  for  the  components  of  strain  at  a  point, 
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du  dv  v 

rd0  dr  ”  r 

7rz 

du  dw 

dz  dr 

dv  dw 

70z 

V  W  J 
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one  obtains,  with  the  aid  of  eqs.  (2.  i)  and  (2.  2),  the  relationship 

le)  =  [q]  la} 


(2.  4) 

(1.  5) 


(2.  6) 


(2.  7) 


where  the  elements  of  [q]  are  listed  in  Table  2.  2. 

The  stress  components  are  related  to  the  strain  components  by  the 
elasticity  matrix  [D],  that  is 

Id)  =  [D]  le}  (2.8) 


where,  for  the  isotropic  case 

*The  symbol  l  }  denotes  column  matrices  while  [  ]  denotes  all  other  matrices. 
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(1  -  y)E _  vE _  VB -  0  0  0 

(1  +  v)(l  -  2v)  (l  +  v)(l-2v)  (1  +  v)  (1  -  2v) 

&  ~  v)B -  - - yi- — 0  0  0 

(1  +  v)(l  -  2v)  (1  +  v)(l  -  2v) 


Symmetric 


E  =  modulus  of  elasticity 

v  =  Poisson's  ratio 

J  ,  E 

G  =  shear  modulus  =  2  (i  +  v  j 


(1  +v)(l  -  2v) 


0  0  0 


GOO 


G  0 


(o)  =  Cor  a9  CJ2  Tr8  Trz  t92)  (see  Fig.  i.  S) 

The  formulation  of  [D]  for  the  case  where  the  material  is  anisotropic 
requires  added  consideration.  Let  1,  2.  3  be  the  axes  of  anistropy  (assumed 
mutually  perpendicular)  and  ay  «3  the  angles  which  define  the  orientation 
of  these  axes.  To  understand  more  clearly  the  significance  of  the  parameters 
a  ,a  ,a  ,  a  step-by-step  rotation  of  the  axes  1,  2,  3  to  their  actual  positions 
is  illustrated  in  Fig.  2.  4.  Let  axes  2,  3  lie  in  the  plane  ABCD.  In  Fig.  2.  4a, 
axis  3  and  plane  ABCD  are  initially  positioned  parallel  to  the  Z-axis  of  the 
element,  ctj  being  the  angle  which  plane  ABCD  makes  with  the  radial  plane 
passing  through  the  central  point  of  the  element  In  Fig.  2.  4b,  plane  ABCD 
is  rotated  an  angle  «2  about  axis  2  to  A’B'C'D'.  Figure  2.  4c  shows  the  actual 
orientation  of  axes  1,  2,  3  arrived  at  by  rotating  axes  2,  3  an  angle  o3  in  their 
own  plane  (A'B’C'D’).  It  should  be  noted  that  if  material  properties  do  not 
vary  in  the  plane  of  2,  3,  then  a3  can  be  arbitrarily  set  to  any  valu*.  say  zero. 

The  stress-strain  relations  for  general  anisotropy  are  given  in  the 


theory  of  elasticity  as 
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(2.9) 


where  the  subscripts  refer  to  the  directions  of  anisotropy  defined  in  the  pre¬ 
ceding  paragraph.  Solving  for  the  stresses  and  writing  the  resulting  equations 
in  matrix  form,  one  obtains 


(o')  =  [d]  t € ' ) 

where 


and 


(o') 
(  € ' ) 


(01  °2  °3  °12  °13  'W 
(£1  e2  *3  T12  ’’'13  T23* 


(2.  10) 
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Table  2.1  The  Matrix  fc] 
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Table  2.  2  The  Matrix  fq] 
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FIGURE  2.  4  Step-by-step  Orientation  of  Axes  of  Anisotropy 


FIGURE  2.  5  Positive  Direction 
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An  equation  connect  the  stress  vectors  to)  ana  »•  can  b  oWabted 

from  consideration  of  equilibrium  of  the  tetrahedron  shown 

for  example,  the  stresses  v*  end  r02  act  on  the  inclined  pUne  ABC. 

Furthermore,  let  (tn>  tJ2,  t23).  <*n'  *12’  ‘13*'  a  31  32  33 
direction  cosines  of  these  stresses  relative  to  the  axes  1,  2. 3.  It  can 
" that  the  components  of  stress  actlnq  on  the  plane  ABC  and  Parallel  to 

the  coordinate  axes  1,  2,  3  are 


R1  *  °lt21  +  T12t22  +  T13t23 
*2  =  T12l21  +  °2*22  +  T23l23 

R3  *  T13t21  +  T23t22+°3t23 


(2.  11) 


.  _  and  t  can  then  be  calculated  to  be 

The  stress  components  Og,  Tf0  ana  tQz 


°9  =  Rl*21  +  V22  +  R3l23 
Tr9=  Rltll  +  R2t12  +  R3l13 
T0z=  Rlt31+R2t32  +  R3t33 


(2.  12) 


,  .  for  0  0  and  t  resulting  in  the  general 

Similar  expressions  can  be  derived  for  0^  oz  a  rz 

relationship  i2.  13) 

(o)  =  IT)  (o') 

ln  which  tx)  IS  a  Stress  transformation  matrix  (Table  2. 3).  1.  can  be  shown 

that  <v  (9. 14) 


The  table  of  direction  cosines  is  given  below: 


where 
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cos  a2  sin(a1  -  0g9') 

=  -  cos  a2  cos(a^  -  0g9') 

=  -sina2 

=  cos  cos(a^  -  0g0')  +  sin  *3  s*n  a2  s*n(ai  ‘ 
=  cos  a3  sintaj  -  0g0')  -  sin  a3  sin  a2  cos(a1  ■ 

=  sin  a3  cos  a2 

=  cos  «3  sin  «2  sin(a1  -  0g0')  -  sin  «3  cosCctj 
=  -  cos  a3  sin  a2  cosfaj  -  0g0')  “  sin  a3  sin(a 
=  cos  a3  cos  a2 
stiffness  matrix  of  an  element  is 
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0  0') 
s 


0  0') 
s 


0  0') 
s 


-  e  9*) 
s 


(2.  15) 


(2.  16) 

i 


in  which 


dV  =  rr  0  z  dr'dO'dz' 
s  s  s 

The  integral  portion  of  eq.  (2.  16)  is  evaluated  by  means  of  the  Gaussian 
quadrature  formula. 


2.  3  Formulation  of  Essential  Matrices  for  Two-Dimensional  Elements 


The  same  relationships  as  in  the  preceding  section,  but  with  terms 
involving  z  and  w  and  components  of  stress  and  strain  in  the  direction  of 
Z-axis  eliminated,  are  used  as  bases  to  derive  the  matrices  for  two-dimensional 
elements.  Th'>  displacement  expression  becomes  simply 


1  r’  0’ 
0  0  0 


r'01 


0 


0  0  0  0 
1  r'  0'  r'0 


(2.  17) 


The  matrix  [c]  reduces  to 
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Table  ! 1 . 3  The  Transformation  Matrix  f TJ 
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(a)  Local  coordinates 
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(b)  Global  coordinates 


TG?JRE  2.  7  Nodal  Coordinates  of  Three-Dimensional  Element 
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(2.  18) 


(2.  19) 


The  formulas  for  components  of  strain  are 


now 


from  which  the  matrix  [q]  is  obtained  to  be 


The  elasticity  matrix  for  plane  stress  is 


(2.  20) 


(2.  21) 


(2.  22) 


Similarly,  for  anisotropic  case 
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CD']'1  = 


'2 
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Sym. 
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*12 


(2.  23) 


With  (o)  equal  to  e<lual  to  toi°2Ti2)  and  °1  deflned 

as  in  section  2.  2,  the  stress  transformation  matrix  is  calculated  to  be 
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,  2 
sin  a 
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cos  a 


-sinocosa 
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cos  a 


sin 

sinacosa 


2  sinacosa 

-2  sinacosa 

2  2 
sin  a-cos  a 


(2.  24) 


in  which  a  =  a.  -  0  0'* 
is 

In  the  element  stiffness  expression  [eq.  (2.  22)] , 

dV  =  rtr  0  dr’dS  (2*  2J 

s  s 

in  which  t  is  the  thickness  of  the  disc. 

It  should  be  noted  that  all  the  matrices,  except  [D] ,  formed  in  this 
section  could  also  have  been  obtained  directly  from  the  three-dimensional 
matrices  of  the  preceding  section  by  simply  deleting  appropriate  rows  and 
columns  of  the  latter  and  setting  values  of  certain  parameters  to  zero.  In 
particular,  to  obtain  [c],  rows  3.  6.  9  and  12  to  24  and  columns  5  to  8  and  13 
to  24  of  Table  2.  1  ere  deleted;  to  obtain  [q]  rows  3,  5  and  6  and  columns  5 
to  8  and  13  to  24  of  Table  2.  2  are  deleted;  to  obtain  [T]  rows  and  columns 
3,  5  and  6  of  Table  2.  3  are  deleted  and  a 2  and  a^  are  set  to  zero.  Thus,  by 
adding  a  few  IF  statements  to  and  generalizing  some  indices  of  a  three-dimen¬ 
sional  computer  program,  the  program  can  be  made  to  work  for  two-dimensional 

cases  as  well. 


2.  4  Failure  Criteria 

The  analysis  described  in  this  chapter  will  be  based  on  any  one  or 
combination  of  the  following  failure  criteria; 


39 


1.  Maximum  principal  stress  theory — 

Under  this  criterion,  an  element  will  be  considered  failed  if  one  of 
the  principal  stresses  equals  or  exceeds  the  strength  (elastic  limit  or  yield 
point)  of  the  material  making  up  the  element.  If  tensile  failure  occurs,  the 
modulus  of  elasticity  across  the  crack  is  reduced  to  zero  and  the  stiffness  of 
the  element  is  revised  accordingly.  The  element  can  therefore  no  longer 
resist  tensile  stress  (no  stress  reversal  is  anticipated)  normal  to  the  crack 
but  can  still  take  tension  or  compression  in  the  other  directions. 

If  failure  occurs  in  compression  (crushing),  an  approximation  is  made. 

To  account  for  the  fact  that  rock  exhibits  ductility  when  subjected  to  high 
confining  pressure,  only  a  portion  of  the  stiffness  of  the  failed  element  will 
be  removed.  A  more  realistic  approach  would  be  to  undertake  inelastic 
analysis  to  obtain  the  actual  stiffness  of  the  failed  element.  This  approach, 
however,  presupposes  the  availability  of  a  stress-strain  curve  which  might 
not  be  feasible  in  most  rocks.  Indeed,  such  curve  can  only  be  obtained 
through  a  triaxial  trest  and  its  shape  will  vary  widely  depending  on,  among 
other  factors,  the  rock  type,  shape  of  specimen,  and  intensity  of  confining 
pressure.  Stress-strain  curves  of  few  rocks  are  available  (1)  but  only  at 
certain  confining  pressure  intensities. 

2.  Maximum  shearing  stress  theory— 

This  criterion  states  that  failure  occurs  when  the  maximum  shearing 
stress  in  a  material  equals  or  exceeds  the  critical  shearing  stress.  The 
Coulomb- Navier  (1)  version  will  be  used.  This  states  that  shearing  failure 
will  occur  when 

-  034i  +  I)  £  2  (2.  26) 

in  which  t  is  the  critical  shearing  stress  found  to  be  between  2%  and  15% 
cr 

of  uniaxial  compressive  strength;  jx  is  the  coefficient  of  internal  friction; 
and  o.  and  o.  are  the  major  and  minor  principal  stresses.  The  above  expres- 

X  O 

sion  takes  into  account  frictional  resistance  to  sliding  along  the  failure  plane. 
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It  should  be  noted  that  If  p  =  0,  that  Is,  there  Is  no  frictional  resistance, 
the  above  inequality  reduces  to 

„  (2. 2 

°1  "  °3  —  ^  Tcr 

in  which  the  left  side  Is  simply  the  expression  for  twice  the  maximum  shear- 
ing  stress  at  a  point. 

Again,  post-failure  behavior  of  elements  will  be  approximated. 

3.  Maximum  principal  strain  theory— 

This  could  be  a  better  alternative  to  the  first  criterion  in  the  sense 
that  here  effects  of  the  other  principal  stresses  acting  normal  to  the  direction 
being  Investigated  are  considered.  To  illustrate  this  point,  consider  rig.  2.  8 
Let  0  be  the  critical  stress  and  ee  the  corresponding  strain.  In  the  first 
figure,  both  the  first  and  the  present  criteria  will  predict  the  same  failure 
stress,  namely,  o.  =  o..  In  the  second  figure,  the  first  criterion  will  predict 
the  failure  condition  o  =  o  .  but  the  present  criterion  will  predict  a  value 
creator  than  oe  since  failure  occurs  when  (Oj  -  Similarly,  a  value 

less  than  o  will  be  predicted  if  o2  is  in  tension. 

Failed  elements  will  be  treated  the  same  way  as  in  the  first  criterion. 


2.  5  Concluding  Remarks 


Derived  in  this  chapter  are  the  different  matrices  needed  in  the 
development  of  the  finite  element  program  proposed  in  this  report.  Discussion 
of  such  other  integral  parts  of  the  program  as  (a)  building  up  of  the  global 
stiffness  matrix,  (b)  determination  of  nodal  displacements  and  stresses,  and 
(c)  calculation  of  principal  stresses  is  omitted  because  they  are  discussed 
in  detail  elsewhere  (2,  3,  5,  6). 

An  equation  solver  proposed  by  Jensen  and  Parks  (6)  will  be  used  in 
the  program.  The  solver  contains  an  optimal  nodal  renumbering  scheme  to 
conserve  sparseness  of  the  stiffness  matrix.  Only  nonzero  terms  of  the 
matrix  are  stored  and  processed. 
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CHAPTER  3 

TWO-DIMENSIONAL  PROGRAM 


3.  1  Introduction 

The  program  described  herein  was  developed  as  an  alternative  to  the 
two-dimensional  program  submitted  as  part  of  the  first  annual  progress  report 
(1).  It  contains  the  following  features  not  found  in  the  earlier  program: 

1.  Element  and  node  information  are  generated  automatically  thus 
avoiding  the  preparation  of  voluminous  deck  of  input  cards. 

2.  Element  sizes  may  vary,  allowing  for  finer  mesh  in  regions  where 
stresses  are  large  and  coarser  mesh  elsewhere,  for  a  more  efficient  solution. 

3.  More  than  one  element  may  fail  during  each  loading  cycle. 

No  change  is  made  in  the  failure  criteria. 

The  element  used  in  the  present  program  is  shown  in  Fig.  2.  1  and 
described  in  section  2.  3.  The  element  belongs  to  the  so-called  "isopara¬ 
metric"  group. 

3.  2  Description  of  Program 

The  flow  of  the  program  is  illustrated  in  Fig.  3.  2.  The  total  program 
is  made  up  of  the  main  program  and  four  subroutines  (complete  listing  is 
shown  in  Appendix  A).  Each  performs  the  following  tasks: 

Main  Program 

—  generates  element  ind  node  numbers  in  the  sequence  shown  in 

Fig.  3. 1; 

—  generates  the  nodal  coordinates,  given  the  radial  coordinates  of  the 

rings  and  the  circumferential  coordinates  of  the  radial  lines; 

—  forms  the  global  stiffness  matrix  in  blocks; 

—  calculates  the  stresses  and  load  factors  (ratio  of  allowable  stress 

to  principal  stress)  of  each  element; 
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-  calculates  critical  loads,  displacements  and  strains;  and 

-  removes  90%  of  stiffness  matrix  of  failed  elements  from  global 
stiffness  matrix. 

Subroutine  DISPL 

--  evaluates  the  unknown  nodal  displacements  by  the  Gaussian 
elimination  method  and  back-substitution. 

Subroutine  MATB 

""  forms  the  matrix  [q]  and  the  matrix  product  [Dj  [q]. 

Subroutine  INTEG 

-  forms  the  stiffness  matrix  of  the  elements,  employing  Gaussian 
quadrature  formulas  In  place  of  actual  analytical  Integration. 

Subroutine  UNIFRM 

—  generates  elastic  properties  of  elements  of  nonhomogeneous  disc  by 
means  of  the  random  number  routine  RANUN  (2). 

RANUN  generates  arbitrary  random  numbers  assuming  statistically 
uniform  distribution.  The  first  generative  number  of  RANUN  can  be  set  to 
any  number  N  by  a  call  to  RANUNS(N).  Both  RANUN  and  RANUNS  are  Madison 
Academic  Computing  Center  (MACC)  library  routines. 

It  should  be  noted  that  in  numbering  the  nodes,  the  center  point  of 

the  disc  is  multinumbered  (see  Fig.  3.  1).  This  is  done  because  each  element 
has  to  have  four  nodes. 

The  well-known  banded  matrix  technique  in  solving  simultaneous 
linear  equations  is  employed  in  the  program. 

3.  3  Test  Problems 

Several  problems  were  run  to  test  the  correctness  of  the  program.  The 
results  of  three  are  presented  here. 

A.  Test  Problem  No.  1 

A  disc,  20  inches  in  diameter  and  1  inch  thick,  is  analyzed.  It  is 
required  to  determine  the  stress  distribution  along  the  line  of  the  load. 
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The  disc  is  assumed  homogeneous  and  isotropic. 

Because  of  symmetry,  only  a  quarter  of  the  disc  Is  considered  In  the 

analysis.  The  quarter  disc  Is  divided  Into  25  rings  and  18  equal  slices  for  a 
total  of  408  elements  and  450  nodes.  The  results  are  shown  In  Fig.  3.  3. 

B.  Test  Problem  No.  2  (Homogeneous  Case) 

A  homogeneous  disc,  3  inches  in  diameter,  1  inch  thick,  and  possess¬ 
ing  xhe  following  elastic  properties 
Modulus  of  elasticity 
Poisson's  ratio 
Allowable  compression 
Allowable  tension 

is  analyzed.  The  disc  is  divided  into  the  mesh  shown  in  Fig.  3.4.  A  listing 

of  the  required  input  data  is  printed  on  page  50.  Some  of  the  results  of  the 
analysis  are  indicated  in  the  following  pages. 

C.  Test  Problem  No.3  (Nonhomogeneous  Case) 

A  nonhomogeneous  disc,  3  inches  in  diameter  and  1  inch  thick,  is 

analyzed.  The  elastic  properties  of  the  elements  vary  as  follows: 

Modulus  of  elasticity  :  5.  5  x  10  to  6.  5  x  10  psi 

Poisson's  ratio  :  0<  23  to  0.  27 

Allowable  compression:  22,  000  to  32,  000  psi 

(allowable  tension  =  5%  of  allowable  compression) 

Double  symmetry  is  assumed  to  avoid  having  to  analyze  the  entire  disc 
which  would  require  considerably  higher  computer  expense.  The  disc  is 
divided  into  the  same  mesh  shown  in  Fig.  3.  4.  Several  runs  were  made.  Dur¬ 
ing  each  run,  a  different  starting  point  for  the  random  number  generator  was 
specified.  The  necessary  input  is  listed  on  page  5 1 .  Results  of  the  analysis 
are  indicated  in  the  following  pages. 


0.  25 

27,  000  psi 

5%  of  allowable  compression 


L 
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(a)  Partial  Disc 

I 


,  {13 

(b)  Full  Disc 


FIGURE  3. 1  Standard  Scheme  for  Numbering  Nodes  and  Elements 
of  Two-Dimensional  Disc 
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I  Start 

_ ^r1 

Number  elements  and 
nodes  according  to 
standard  scheme. 

=Z  * 


Carry  out  elastic 
analysis  due  to 
1000  lb.  diametric  load. 


r  ♦ 

Calculate  load  factor  (LF) 

of  each  elemeat.  Let  LF 
denote  critical  LF. 

^  T 

Consider  elements  with  LF 
within  X%  of  LFcr  as  failed. 

X  is  part  of  input  data 

_ r  1 

Calculate  critical  load  (1000.  x  LFcr) 
and  corresponding  vertical  displace¬ 
ment  at  point  of  load,  horizontal 
strain  at  center,  etc. 


No 


—i 


Deduct  90%  of  stiff¬ 
nesses  of  failed 
elements  from  global 
stiffness. 


Note:  Load  factor  is  the  ratio  of  allowable 
stress  to  principal  stress. 


FIGURE  3.  2  FLOW  CHART  OF  TWO-DIMENSIONAL  PROGRAM 


0 

y 

u 


FIGURE  3.  3  Stress  Distribution  Along  Line  of  Diametral  Load  (Test  Problem  No.  1) 


INPUT  DATA  FOR  TEST  PROBLEM  NO.  2  (HOMOGENEOUS  CAS 


INPUT  DATA  FOR  TEST  PROBLEM  NO.  3  (NONHOMOGENEOUS  CASE 


NOTE:  In  this  set  of  data,  starting  point  of  random  number  generator  is  N  =  3962185381  (Card  No.  5) 
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3.  4  Discussion  of  Results 

The  results  from  Test  Problem  No.  1  clearly  demonstrate  the  theoretical 
correctness  of  the  finite  element  solution  formulated  in  this  report.  As  can 
be  seen  from  Fig.  3.  3,  the  elasticity  and  finite  element  solutions  differ  by 
only  6%  for  and  2%  for  o  at  the  center  of  the  disc.  By  solving  the  same 
problem  several  times  progressively  decreasing  element  size  each  time,  the 
finite  element  solution  was  also  found  to  be  convergent. 

The  discs  analyzed  in  Test  Problems  No.  2  and  No.  3  are  models  for 
the  Valders  limestones  studied  extensively  in  the  experimental  phase  of  this 
research  (Chapter  1).  The  elastic  properties  assigned  to  the  nonhomogeneous 
disc  are  actual  range  of  values  obtained  from  experiment  .while  those  assigned 
to  the  homogeneous  disc  are  the  mean  of  these  values.  The  failure  pattern 
predicted  by  the  finite  element  solution  (Figs.  3.  5  and  3.  6)  agrees  favorably 
with  the  actual  failure  patterns  of  Brazilian  tests.  The  critical  load,  however, 
appears  to  be  underestimated — 5730  lbs.  and  4400  lbs.  for  homogeneous  and 
nonhomogeneous  cases  compared  to  the  6940  lbs.  average  obtained  from  the 
experiment.  This  seems  to  indicate  that  either  the  90%  factor  used  to  deduct 
the  stiffness  matrix  of  failed  elements  is  too  high  or  the  failure  mechanism 
assumed  in  the  program  is  not  altogether  realistic  or  both. 

The  load-displacement  curves  appear  to  follow  the  same  general  shape 
regardless  of  material  characteristics  cr  starting  point  of  random  number 
generator  (see  Figs.  3.  7  and  3.  8).  These  curves  compare  favorable  with 
those  obtained  from  the  experiment  at  the  early  stages  of  loading. 

3.  5  Conclusions 

Based  on  the  several  computer  runs  made,  the  following  conclusions 
are  drawn: 

1.  Even  in  its  oversimplified  form,  the  finite  element  solution  is 
capable  of  predicting  the  actual  failure  pattern  of  Brazilian  tests. 

2.  With  a  few  improvements  in  the  failure  criteria,  there  is  a  strong 
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possibility  that  the  critical  load  and  the  load-displacement  curve  can  be 
accurately  predicted  as  well. 

3.  In  nonhomogeneous  cases,  the  shape  or  the  load-displacement 
curve  is  not  affected  significantly  by  the  choice  of  the  starting  point  of  the 


random  number  generator. 


Figure  3.6  Progression  of  Failure  in  Test  Problem  No.  3 


5000 


FIGURE  3. 
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APPENDIX  A 

two-dimensional  program  listing 

The  computer  program  presented  here  Is  written  In  Fortran  V  language 
specifically  for  the  Univac  1108  computer  at  the  University  of  Wisconsin  - 
Madison.  The  program  Is  capable  of  handling  all  sor*s  of  static  two-dimen¬ 
sional  loadings,  not  Just  the  diametric  loading  of  the  Brazilian  test.  Zero 
or  nonzero  displacements  may  be  assigned  to  any  node  making  possible  solu- 

tions  involving  portions  of  the  disc  only. 

The  program  requires  usage  of  three  auxiliary  storage  tapes  or  drums 

designated  by  the  numbers  10,  11  and  12. 

The  listing  is  complete  except  for  the  random  number  routines 
RANUN(R)  and  RANUNS(N).  These  subprograms  are  provided  by  the  Madison 


Academic  Computing  Center  (MACC). 

A  description  of  the  input  and  output  parameters  required  in  the  pro- 

gram  is  given  below. 

InDut  Notations: 

1.  (a)  NCHO 

-  Printout  code;  0  if  only  stresses  at  failed  elements 
are  to  be  printed  out;  nonzero  if  stresses  at  all 
elements  are  to  be  printed  out. 

(b) 

NWCL 

-  Node  at  which  radial  displacement  is  to  be  computed. 

(c) 

NST1,  NST2 

-  Nodes  between  which  strain  is  to  be  computed. 

2.  (a) 

EM1.EM2 

-  Range  of  values  of  e’astic  moduli. 

(b) 

CS1,  CS2 

-  Range  of  values  of  allowable  compressive  stresses. 

(c) 

PR1,  PR2 

-  Range  of  values  of  Poisson's  ratios. 

4 

(d) 

CTRAT 

# 

-  Ratio  of  allowable  tension  to  allowable  compression. 

(e) 

T 

-  Thickness  of  disc. 

(f) 

TRN 

-  Maximum  percentage  variation  from  largest  load  factor 
for  element  to  be  considered  failed;  0  if  only  one 
element  is  allowed  to  fail  each  time. 
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(g)  DIA 

3.  (a)  IRN 

(b)  NIK  fl) 

4.  (a)  NELEM 

(b)  NTNN 

(c)  NNSD 

(d)  NNCL 

(e)  NE 

(f)  NS 

(9)  NC 

(h)  NER 

(i)  NCOD 

5.  N 

6.  RADN(I) 

7.  TTAN(I) 


-  Diameter  of  disc. 

-  Total  number  of  elements  in  vicinity  of  applied 
loads. 

Number  of  Ith  element  in  vicinity  of  applied  loads; 
specified  if  run  is  to  be  terminated  upon  failure  of 
element;  0  if  run  is  to  continue. 

“  Total  number  of  elements. 

*  Total  number  of  nodes. 

■  Total  number  of  nodes  with  prescribed  displacements. 

Total  number  of  nodes  at  which  concentrated  loads 
are  applied. 


-  Total  number  of  slices  portion  of  disc  involved  in 
analysis  is  divided  into. 

Total  number  of  equal  slices  whole  disc  is  divided 
into;  0  if  disc  is  not  divided  into  equal  slices  and 
only  a  portion  of  whole  disc  is  Involved  in  analysis; 
3fly_ negative  number  if  disc  is  not  divided  into 
equal  sttces  and  whole  disc  is  involved  in  analysis. 

Total  number  of  rings  plus  one  disc  is  divided  into. 

Least  total  number  of  similar  consecutively  numbered 
elements;  1  in  nonhomogeneous  and  anisotropic 
problems.  This  parameter  allows  the  generation  of 
stiffness  matrix  of  similar  elements  only  once. 

-  Total  number  of  load  cycles. 

First  generative  number  of  random  number  routine. 
(Note  that  this  information  is  supplied  only  in 
nonhomogeneous  problems. ) 

-  Radial  coordinates  of  rings  starting  with  0. 

Circumferential  coordinates  in  degrees  of  radial  lines 
dividing  disc  into  slices.  (Note  that  this  informa¬ 
tion  is  supplied  only  when  disc  is  divided  into 
unequal  slices. ) 
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8.  (a)  MR 

(b)  NODC(MR) 

(c)  C  LOAD  (MR,  I)  - 

9.  (a)  NODBOX) 

(b)  DISPOX,  I) 

(c)  JJ 

(d)  n 

Output  Notations: 

RADIAL 

CIRCUM 

SHEAR 

PRNCP1,  PRNCP2 
MAXSHR 
LD  FAC 

All  other  output  notations 


1  for  1st  card;  2  for  2nd  card;  3  for  3rd  card,  etc. 

Number  of  MRth  node  with  externally  applied  con¬ 
centrated  load. 

Magnitude  of  concentrated  load.  (I  =  1  in  radial 
direction;  I  =  2  in  circumferential  direction. ) 

Number  of  IXth  node  with  prescribed  displacement 

Value  of  prescribed  displacement;  200.  if  displace¬ 
ment  is  not  prescribed.  (I  =  1  in  radial  direction; 

I  =  2  in  circumferential  direction. ) 

Total  number  of  nodes  which  exactly  have  the  same 
prescribed  displacements  as  NODB  and  which  form 
an  arithmetic  progression  with  NODB;  0  if  none 
such  nodes. 

Interval  of  arithmetic  progression. 


Direct  stress  in  radial  direction. 

Direct  stress  in  circumferential  direction. 

Shearing  stress. 

Principal  stresses. 

Maximum  shearing  stress. 

Load  factor  (ratio  of  allowable  stress  to  actual  stress), 
are  self-explanatory. 
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f*Int  iuii7.rai.cp7 

rm?  FORK  AT  »  1 3HIIFC  rSFUO’S  RATIO  :  .F|u.*.4H  TO  .Fju.rj 

P .? x n r  Iona .r si .r S5 

1MUG  f OtOK AT  I  71H)irc  TT  TCAL  CbMFRFCSTVr  STRFSS  =  . 
nrf :u,2  .dm  to  » ?Ft im. i  « 
r«TNT  UHir.  fTRAT 

1  'T'*0  f  0 .0  K  A  T  f  77HI1C.7T  TTCAL  TFNSIGN/r.OMrRFSSION  RAT  TO  =,F*.J1 

PkTf'T  lllllul 

ini»  fohmat  « ? »i H* i r h r r k s %  of  orsc  =.F5.?i 

FOrUT  n  il.  OTA 

1*511  FORMAT  *  1  r»H»:OTAMF  tcr  OF  Olsr  :,Fc,£) 
rfat-  ic<ni.  :*5.n.  cntx  *  n.  r  =  i»  trm 

nrA.O  lo.ll  « NrL  F  •!  «'i  TNR  »  NNSO  .  NNfL  .  OF  .  N  %.  .  OF  R.  NC  OD 

1M'1  roKKAT  *niri 
NC  1  =  *.'C  - 1 

PR  TNT  If  1  Z  .  \'F l  rM 

I'll’  r  O.IMA  f  *7?,Hi  TOTAL  OF  F|_»rw-^rc  -  #tt| 

FKTKT  l*i !  .*  •  NTf.’fv 

i  *•  i  *  format  *  ? .  ho  total  Njui.a  of  no.ofs  =  »e*i 
trim  iiin.f;f.:,c: 

I'il'i  FORMAT  (  *i  1  U',  a  T  jC  )•>  fokttom  of  o  isc  rs  GTwroFo  into  . 
UTCfS  AKF  *  I'.IH  RT^CSI 
Tr  t  A >  *  F  M 1  -  r »!  <  .o:  113*. 173.1  21 

i  •  i  «*  r  ?  r  i : r « .  u 
117?  FORMAT  ITMI 

TALL  RAMJOS  *M 
00  1  Ml.  :=i.%flfs 
T .'Ll  U •*.  T  frm  t  r  .M  l  •  f" e  •  f  »•  S  I 
FALL  bJTF  RM  »?PI  .?■*<  .PRSI 
FALL  U  -J  r  F  (  f7  i  .  Cf  t  •  C7  r. ) 
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«■  r . 
r  ■». 

T  F  . 
e9. 

I  n. 
r  t . 

I  5  . 
f  ’ . 

'  4  . 
Fc  . 
f  r . 
r?. 
e  P. 

n. 

7f  . 
7  l  . 
75  . 
•»’. 
7b  . 
7r. 
7f  . 
7’. 
7  r  . 
7  2. 
fill. 

ai . 

p' . 
a**. 

P4 . 
PS. 
£f  . 

a?. 

5f  . 
fi9, 
2t'. 
91  o 
r: . 

27. 
9». 
9r-. 
°C  . 
9’. 
op. 

9". 
I  ItM  . 
lit!  . 
1  <■'  . 
Itl7. 
1  HU  . 
lur. . 
1  lif  . 
1U?  . 
I  i*P . 
!!!'•. 
1  1  •«. 
111. 
115  . 


r  re’. :» ms 

F-|(  7.2  l  =  POr, 
li  t.  PPf!.  Iir-CS^ 

PRINT  l  If  **  .FO  ('  •  1  I  ."D  c  ?  1  »!  » 
lirn  TOM  AT  <  IFf  l  5  .f  ,  IF  t  12  .f  I 

nn  to  i  if 

1**  ro  If*  Tri.NPLP« 

P9C T • 1 l-Frtl 

rner. : » =r  k i 
ini  roc  :.i  »  =  *csi 
105  ri;r\'T  F997 

roo?  ro.VlAT  1  3H0FLf,.v*r‘*  T»S  »  •  1?MC  0SNP3  N0DFS.3k.3H  FOFMFMT. 

If  *  .  l’HCOFNFO  N  0  D  F  S  t 
701,:  FOR  “A  I  < e  TF  I 

rr  eNF-wsi  9t  «  oi.  l  *r 
12F  TF  C  MS  191  •  «')•  915 
«s».  00  ? ii I.»  Trl.NFLFM 
NCCt  T. 1 »  =  I«  < I  -1  l/NF 

None :. 21  =nooc  r .  1 1  ♦  i 
\r,ce  i*3 » = Mor> « ra 2 1« !  «ne 

M1D(  1.41  =N0O(  T.  *1-1 

5iim  fonttnuf 

N '»»  :  CNF*  2 »  •  2 
oc  to  i  «;i 

f*l  no  192  lri.Nn.FM 
None  t,i  i  =  i 
None  Tf  21 r T«  1 
Nooe  t.  !  i-Noae  :*2  mnf 
NODC  T. 41 =  NUDC  I. !l-l 
199  C  ON TTNUF 
NPlI  =N  £•  4 

OC  198  i  =  nf»nflf,v  •?;? 

Non  r.n  =hodi  r.  1 1  -nf 
nod e  it 3  i  =  Noner.’»  j-:jf 

199  CONTTNUF 
l *»:  ngznflfm/z 

N  C 1  5  C  N  fl  f  N  ♦  II/? 

no  73  r=i tSG 

Krl  ♦NCI 

print  tf  1 1 1  e None  r » .j  i •  j- i  t«i » tK* e nooiki  j i  t  j=i  »«i  1 

?'.  COMTNUF 

■♦f  FORMAT  f*  If.  16  I 
IF  e  NO-NCI  I  *9t  1  9  7i  v:i 
»«;  PRINT  7nlufNGl .  INODCNGl  »  Jl  t  J-l  *4  I 
f  *  •  •  *  •  G  ? N  rS  AT  F  NOPAL  COO'  n  IN  a  T  r  s*  ••••••••••••••••••••*••• 1 

197  Or. AT,  1UH2  f  <  9A0:.e  T  1. 1=1  fNC  I 
1  It'll  r  oRN  AT  C9Pf.  .21 

rr  e nf_ -ns » no  .  in: .  j  7? 

1*7  Ic  C  NS  t  i  f:  2 »  l  ?4 »  l  94 
12  7  NF=NF-l 

12b  no  l^r.  ;  = : .  ntnn 

NDUM= i * «  T- 1  » /CNF* 1 > 

r.  tpe  r  j  =sAr:e:er.nuNi 

19T  CON TTNUF 

f '  F  :*  F  ♦  i 

rr  e  ns'2*1  i  #:m  •  5mz 

-til  pr.r  i  '.,i,  r  •  f  T  T  AN  (  I )  •  I  =  1 1 N  F  • 
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1 1  r . 

114. 

1 1  p. 
110. 
117. 

—  1*  — 

in. 

?l'  * 

1 1 «. 
l  ?*»  • 

12  1. 
1"*? . 

2n4 

12  7. 

r: 

12“  . 

12  E. 

nr 

1?f. 

12  7. 

1UI  r 

in . 

1  2  c  • 

llliC 

171*. 

nt. 

r*  *  *  * 

1  ?7  . 

2111m 

i 

1  ’4  . 

iwi: 

111. 

c  *  *  *  • 

1  7  F  . 

rnni 

PI. 

m. 

:uii4 

no. 

1  4(1. 

m. 

IM'14 

14". 

-  r  •* 

i. 

14’. 

144. 

141. 

1  4F  . 
14"’. 

*14 

1  4  A, 

:r  r 

14  1. 

1“'. 

2r  o 

Ill  . 

-no: 

:r* . 
ir  T. 

n». 

'Cf 

A  ■  • 

ir< . 

r*  ♦ 

1*7. 

l"’ 

1  c.  0  . 

r ..... 

1*  r>. 

r 

1  ‘III . 

c 

u  1 . 

IF*  . 

If  7. 

ns. 

1  E  r  • 

r 

Kr . 
1C7. 

IFF. 

IF  F. 

1 M4 

r.o  to  :nt 
*«Ni  =  NS 
nu  K  =:  T  E  tl .  ✓  HN<^ 
nc  "l'3  1“  1  » Nr 
t>I-T 

TTANtrir(AI-l.  I  *0 «JM 

nG  10£  iri.fjr 
7T.z  T«  Nf  •  (  NC  -l  I 
rDU»*rITAN(  n  ♦  .I<1  74 

a  a  102  j=r.rt,».f 
T  T  A(  J I  -T  PuM 
CONTTSUF 
F  W  T  M  T  11, 1  £ »  (VADNiri  #  r  =  l#NC» 

FORMAT  ( 1  3K.SRA0  TAL  COG  RD  INA  TF  %/ (  G  *  *  3F 1 1! ,  2  I  1 
PUMT  1H  16.  II  TAM  lit  I=1.NF  I 

’■jliriAT  <  7  3  Hil  C  T'TCUHFEflt  NTTAL  COOROIMTFS  OFGRFCS) 
1/  <  f  X  .  GFK.  .ZI  » 

•  TliPL'T  MAGN  ITtJf'tS  Qf  f  PS  T  EN  T  f.  A  T  F  Q  LOAOS  If  AM**** 
TF  I  MNCL  I2UU1  *  ‘ l!IU  *  2Ulif> 

RFAT:  ilir.».W??.NPOCCVRI«  (CLGAniMP,  Tl»  I  =  i»  cl 

FGI^AT  ITc.TI.JFIu.II 

Tf  I  MR  -NNCL  »  cl'llfi.  JfiUI#  2  till  1 

•  T‘JPi;T  VALUfC  or  sFFC  TF  rr0  nrSPLACFKCKTSA**-******! 

r  x  ='.i 
r>- tx« i 

a  r  An  i  mm,  no  op  ( rx  i  •  in  r?  f  i  rx  *  n  » r  =  l » 2  i .  j  j»  it 
F0Rmat  CT£*cF£,l,r?*T3l 

C'UM  Z  T> 

TF  (  T> -NNSn  1 2£  2  * 2ltU 2  f  c'UU2 

TF  I  Jfc:l  J  HU  4  »  2 Mil 4  *  2  r  4 

Ho  21G  T  =  1 1 0 J 

T« =  r  •  T I*  NODF ( NOUKI 

T  *z  lx«  1 

f.ftJPt  IX  I 

I'll  2F5  jri,; 

01  SF  I  rx  ,  0»  r  r.lSF  |UruH»  J» 
f  o  in  t  r  fl  U  F 

l’fi5 »  Z  till 2 


r  r 

\t r'  - 


i  rx  -nmpi  n  im* 
mi; 


*m= :  -Xtr 
rC  1 M  T=I.KV 
r(  I  1  =  0. 

0!i  ! !  3  J? 
c  I  r  .  01  rfl. 


ft  irr^rr^  nats-tcf^  tc  ?rc;?*. »»»•»»«»•••*' 


- 1  .Nn 


TRFATTNfi  TUG  FLOCKS  FACh  T  rvF  mmmmmm.m 
SFAprn  rp.T  rLF-.r?,-T  THA  r  CONTRTF  G’T«n  TO  UPPFR  STOCK 
A  NO  C-  r  N  r  F  A  T  r  ITS  attffnfrs  ant  lpac  katptcf<: 
RF  A  r NO  in 
R  C  tl  T  N  ?  1? 

L  f  N  = -l 
t!  FT:  r  r 
MlNrNP*.*  1 
KL  -IT.'Rf:-  11  ♦  1 

N  1=  N^K*M  l 
M2  “N T  «r;i  *» 
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1  7ft. 
171  . 
17?  . 
1 7  "*  • 
17b. 
I7e  • 
1  7F  • 

177. 

178. 
17®. 
l  crs. 
181  . 
i  e?. 

187. 
18b. 
1  8C  • 
186. 
187. 
1  88. 
18®. 

190. 

191. 
18?. 
191. 
19b. 

nr. 

1  9F  . 
197. 

1  96. 
199. 

2  00. 
2»!1  . 
2  0?  • 
?fP. 
2  U*»  . 
20*  • 
2  06  . 
2  f  1 7  . 
2UB. 

209. 

210. 
Ml. 
212. 
217. 

zib. 

21*. 

216. 

217. 

218. 

219. 

220. 

221. 
m  ^ 

c  c  t  . 
227. 
2Tb. 
229  • 
226. 


00  10*  Trl.NFLFM 

rr  inooit.i  mns.ios.iii*; 

1  !if  TN:T 

00  107  J=l.b 

IF  INOOl  I.  J1  *2-NTt  108.  108.  107 
107  CONTTNUF 
r.O  TO  10* 

11.  G  r  13N00II.11 
T2rN0  0l  1.21 

Mr  .**1RAQ1T2  I  ♦  RAO  I  II  »  1 
X  2  r . 6 *  I R  AOt T2I -PA01 Til » 

72  =  .*=•«  TTA«  T2  I-TTAI  T1  I  I 
T  F  1721206. 206.207 
MG  V2=,*»I6*28M9-TTAII1  1  I 
207  rirTT  A(  Tl  M  V  2 
Fr  °D1 1 » 1 1 

TRIP  01  1.21  _ _ 

ff RrFOl T.3  1 

URITF  1 121 IN » X 1.X2.T1.T2.F.PR.CCR 
LTNCrl  l-l  1/NFR 
T F  1LVNC-LVN lllb.il U.llb 
lib  LTN=LTNC 
CALL  INT  rG 

C  A 00  fLFMENT  "TTFFNFSS  MATRIX  TNTO  UPPFR 

C  TRIANGULAR  PORTION  DF  SlT.JI 

110  OC  219  Kri.b 

00  ? 1 3  Jr  1 . ? 

KNrtK-1  »*2«J 

KM  =1  NOOl  t.KI-ll*2«J-INfiN-ll*NL6 
00  212  L=l.b 

TF  fNOOl  T.L1  -NOOl  I.K1  1212.  Tib. 21b 
21b  00  21*  M>r 1 *2 
LNrCL-ll*2*KX 

L Mr  I NOOl T  «L l-NOOf I.K1 M2«MX-J«1 

TF  1LMI219.21T.21 2  .  . 

212  81KM.LMirSlKM.LMMSLlKN.LNI 
21*  CONTTNUF 
217  CONTTNUF 
219  CONTTNUF 

NOOl T.l  ir-NODl T.l  1 

1 1J9_  CONTTNUF  _  _  _ 

C...**rNTROOUCF  CCNCFNTRA  fED  LOAOS  TO  firi»»»MMM».»M»  ......... 

TF  1NNCLI  HI.  10.  2 
9  00  b  trl.NNCL 

TF  1N00C1  IM2-NL  1 b.F.9 
2  TF  1  N00C1 1  1  *2 -NT16.6  .b 

6  00  7  Jr  1, 2 

TF  (CLOAOl T.J  1  18.7.8 
8  KMrlNOOCl  11-11*2  •  J  -  INFN-l ) *NL 6 
FI  KMlrFlKMMCLOADl  I  .J! 

7  CONTTNUF 
b  CONTTNUF 

C»*"»  TNTRO  DUCF  90UNOART_CONOtT  IONS  TO  SIT.JIANO  nriMM.MMM.M 

11)  OO  lb  in.NNSO 

TF  1N00B1  lM2-NLUb.l6.15 
1*  TF  l N008 111*2 -NTI16.16. 17 
16  CO  15  J - 1 , 2 
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5  2  ?• 
2  2  8  • 
22*?. 
2  in. 
221. 
215. 
222. 

214. 
275. 

215. 
2  2  7. 
2ia. 

224. 

240. 

241. 
545. 
2  4  2. 

244. 

245. 

246. 
24  7. 
244. 
24®. 
25fl. 
211  . 
252. 
252. 
254. 

2  55. 
256. 

2  5  7. 
254. 

2  5  9. 
260. 
261. 
262. 
262. 
264. 
26  5. 
266. 
26  7. 
264. 

2  C  9. 
27r . 

2  71. 

275. 

2  72. 
274. 

2  75  . 

276. 

2  77. 
578. 

2  79. 
280  • 
281. 
292. 

2e*. 


rr  c  idd.-absidtspi  t.  ji  i  ue.i8.25 

2  5  KM=!NOOOm-l  I  *2  ♦  J-f  N8N-1  I  *NLB 
FfKMl =DTSPC T.JI 
SC  KM.l 1=1. 

KKrKH 

00  19  K=5.N«M 
KK  =KK  « 1 

P(KKI=FCKK I -SC KM.KI *Ff KM| 

19  SCKH.KlrO. 

KK=  KM 

00  20  _ _ 

KK=  KK-1 

TF  CKKI2i.51.22 

52  F(KKI  =  F(KK  l-SCKK.M*  I  *0  I  SP  (  I «  J  I 
2(1  SCKK.MXI  =0. 

51  C  ON  TlNUf 
18  CONTINUE 
GO  TO  14 

17  KRK=NOOB(  II *2-NT -NBA 
TF  (KRK  121.51.14 
2?  DO  24  J  =  l.  2 

TF  C100.-4R5f0rSPCT.jm54  .54.56 
2  6  KM  =  CN ODBC  II  -  X.l_*5  » J  -  (NBN.-1 1  •  NL  9 
LRL-KM* 1 -NBW-NLB 
TF  CLRLI 58.58.54 
58  KK=  NL8  *  1 

L  IM 1 =KM-NL  B « 1 
90  57  ML=L IM1 .NRW 
KK=KK-J 

FfXKlrFfKK l-SCKK.Ml I'OrSPf  T.JI 
S  CKK.MU  =u. 

57  CONTINUE 
24  CONTTNUF 
14  CONTTNUF 

C*.***?MR.I.TF_  UJPPFR  .BLOCK  ON  TAP  F  4^0  _S  HITT  UP  10«F»  BLOCK **«»•••»••« 
URTTF  C 10  I  CFCNI.CSCN.MI «M=1 «NB  Ml « N=  1 . NLB  ) 

TP  CNTNN*5-NL6I  50.  c(3. 11 
31  TF  f  NTNN*2 -NT  If 0  *50 .15 

12  00  I?  1=1. NLB 
LL=L«NLB 

raj=FLu» 

FCLLUO.  "  . 

DO  32  MLU.N9M 
SCL .ML l  =  SCLL .ML  > 

13  SCLL.MLI =U« 

60  TO  104 

6*  •  •  •  •  COMPU  T  E  AND  PPTNT  NOPAL  D  IS  P  L  ACF  KFNTS  •»»»•.*•»•••••••»**....  . 

51)  CALL  OISPL  - -  —  "  * 

IF  CNLBI 598. 995. 138 
118  DO  40U  MJ= i . NC QU 
FRTNT  2999. MJ 

5939  FORMAT  ClUHOCTCLF  NO. .Til 
CRLFr  lDUflilin  • 

TF  CNCH0n6.95.9S  ~~  . . 

95  PRINT  1000 

lOOn  FORMAT  <  3  8  HU  STRESS  ’S  A  T  ELEMENT  CENTRAL  POINTI 
PRINT  *002 
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1  54. 
2Cr. 

2  TF  . 
29?. 
2  SC. 
280. 
2  9(1. 
20 1. 
2  92  . 
29?. 
2  9*i  • 
2  99. 
✓  ‘'r. 
297. 
2  OP. 
709. 
2((|l. 
llll. 
2(12. 
’ll*. 

2  li*l  • 
309  • 

3  lif  • 


P.i^  fORPAT  (Sri  fL(MfMt7H  k  AO  IA  t *  £*  *  P  HC  TRCU  **«  6*  *  £  H  SHEAR. 
1  H  X  .  f,  HpRNCP  1  .6X.SHFRNCP2.CX  »6HMa  XSHR.6X.GHL0  FAC  ) 

°r  rfxtnd  12 

00  4F5  Ml=l.NFLfM 

OfAO  (  121  TN.X  1.  *2.yi.r2.P.PR.CCR 

TF  (  NGOI  TN.4  I  1465  *92 .92  . 

22  N90(  TN.  llrAPSfNGDdN.  Ill 
TCR  = -CC  R*C TRA  T 
S  CR -- CCB 
00  4F 6  J=1.4 

Kf  =1MNOO<  IN.JI*  2-ll/Nl  B 
L  =  <  NOOI TN.J  I  -1  I  »2 -NLR • I KP-J » 

ITH'=2* J 
LTH1=L  TM7  -1 
ro  An b  ktl  rvi.L  w.i 
L  =  LM 

Al'E  r m  =S(L  »KF  I 
466  CONTTNUF 

00  81  1=1,4 
OLf  r  i=o. 
oicr«m=i'. 

00  flu  K  =  1  .7  .2 
KK=K«  1 


2*17. 

Kl=XX/2 

2  liR. 

OL  (  r  1  =01  (  n  «  C  (  I*K  1  i  *r  IK  I 

300. 

90 

0L(  r«4  l=HL  <r«4  Mf  (  r  fKl  |  •FIKKI 

31(1. 

81 

CONT  TNUF 

311. 

X( 1  1  =  6 A ( *  | 

312  . 

mi  -  G6(  i  i 

313. 

M=1 

314. 

Nil 

319. 

CALL  HA  TR 

31  r. 

oo  ri  r=i,j 

317. 

2(11=0. 

218. 

00  FU  J= 1  •  8 

310. 

RO 

2( r  i=zi r  i  «eo( r t  ji «ol ( j i 

32«. 

F.l 

CONTTNUf 

!?i. 

OJHl=SQRT(.29*(2(l)-2(2ll**2*Z(3l 

22  2  . 

DUM?  =  .;*(Z<  iJ  «? (2»J 

32  3. 

T5=0UM2  «0UH1 

32  4  . 

CS=0UH2-CUH 1 

329. 

S  S=  OlJHl 

22  r  • 

• 

72ITN.il =7( I | 

327. 

221  TN.2  1=2(7  1 

32  R. 

72 (TN. ?| rZf 71 

320. 

721 TN.4  1=  TS 

32(1. 

22 (TN.EI =C9 

331. 

721  TN* 6  l=SS 

33?  . 

7R  =1 . UF-f 

331. 

TF  ( TS-7RI62 .92 .6 3 

2*4. 

C  * 

TF  (  Cf.  *2R  1  E9 •  f  7»  F- 7 

339. 

F5 

TF  (CCR/CS-TCR/TSI87.67.67 

23 f  . 

87 

HLF=f  CR/TS 

337. 

K€F( TN»=i 

2  *  r. . 

CO  TO  £4 

!  30  • 

52 

HIF  =  CCR/C  i 

3411. 

KfCI  "01  =li 

21 
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3  A  1  . 

LA 

3A?  . 

3C 

3  A’  • 

r  e 

2AA  . 

3Ar  • 

69 

JAG  . 

JA  7. 

7A3. 

96 

2  A  5. 

A 11 02 

3911  • 

97 

321. 

f  6 

3?  2  . 

55*. 

2C  A  . 

Af  c 

*••9. 

*rf  • 

1  ?  9 

32  7. 

1 2  f 

39  8. 

1  ?  2 

329. 

3  F  1 1  • 

361. 

3f  5  • 

267. 

A  Oil  2 

3t  A  . 

3E2. 

2FF  . 

367  . 

ire. 

*69. 

3  70. 

271. 

27?  . 

273.  IlSt 

?  7A  . 

27  2  . 

2  7f  . 

277.  76 

2  7*. 


fF  lAf-SCSSI -ZRIC9.F  2.89 
rr  1  AiISCSCR/SSI-HLF  163  .€9»E9 
HLCr  ABSCSCR/SM 
*FF  1  TN  »  =  -1 
rCR  =  lUUP.*W-F 
FLOFJ  TN  »  =  HLr 
TF  IN  CHOI  9  7.  98*9 
PRTNT  AGU2.TN.ZU 
FORM A  T  ( I9»  8F1?.5 
TF  CHLF-C  RLF  I  6E« 

CRLFrHLF 

PCRCrFCR 
NO-  TN 
CUNT  TNUF 

TF  C NCHO » 12c * i ?? • 122 
TRINT  12C.FCRG 
FORMAT  117 Hfl CRITICAL  LOAO  - 
NBKHNHCL-U  *2/NL3*l 

LHi<  =  NUCL*2-l-«NBK-l  l*NLn 

M  0K -NL  B*  L  BK 

nsP=PCRG/io«io.  •sfinM.NftK  i 

FORMATC12MI*  VERT.  >9  T  Sr  L  •  AT  POINT  OF  LOAD  =.FU.5I 

NBK1U  NST 1- 1 1 • ?/Hlp  *1 

N9K2  =  <NST2-1  I  *2/NLP*l 

L6KirNSTl*2-l-INF^l-l»*NLfi 

L9M2-NST2 *2 -l -INP*2 -1  I *NLO 

MBK1=NLB«L&KI 

MBK?  =  NL«*LRK2 

RTRrtStHRK2.NBM2l-SIMBKl.NPKl  I  I /IP  AO  I  NSW  I 
1  -RAOtNSTl  I  I  .PCRG/iOllU. 

PRINT  A1P2.  SIR 

FORMAT  t?OHt:STRATN  AT  CFNTFR  =  .Fl<.,» 

JFNrfl 

r* U  M67  MK r  1  » NF L Fm 

IF  CNOOI MK» A  I  I  AC  7. AC ?» 7f 
TRUNrtFtnFtNOI-FLOFtMK  » l/FLOFINOI 
IF  I  ARSt  TRUNI -TPM  77.  77.  *567 


279. 

?  P 1 1 . 
Irtl  . 
3  BZ  . 
38’. 
*8«i. 
28e  • 
3  8f  . 
387  . 
3  B  8  • 
3  8°  • 
3911. 
391. 
2e7  . 
3*3  ■* . 
3  9A  . 
Znr, 
2**F  . 
297. 


7  7  JFN  =  JFNU 
KFLt JFNI rMK 
DO  AF3  MC-1 • IRK 

TF  f  MK-NTKt  MCI  »  AC  B*  7B»*GB 
7  fi  MAMCHA=MK 
GO  TO  9° 7 
Af.8  CONTINJF 
(f f  7  CGNTTNUF 

TF  t  NCHO  I  9A  .2  2*1  »°A 

AOUl  FORMAT  f'HIJTHF  FOLLOWING  FLf.MFNTS  H»VC  FATLFOI 
CO  TO  122 


FTRFSSES  AT  fatlfo  FlFMFNTSI 


<3A  PRINT  2*1*!  3 
llO»  FORMAT  I  3  A  HI' 

fu.A  FCftMAT^BH  FLFr.£M.7M  R  AO  TAL.  6X  .  6HC  TP  CUM,  C  X  .  EH  SHFAP. 

lCX.6HFRNCFl.5X.CH?RNCF2.6X.6HHAXSHR.6x.l2HFAILU«r  "0  IF  I 

172  no  AF.  S  Mn  =  l,  JFN 
NOrNFLCM.il 


69 


too. 

«i  on. 
MIU. 
Ml)'. 
MU?  • 
MOM  . 
MUf  . 
Mijr. . 

Ml. 7. 
Mllfl  . 
Mil**. 

mn. 

Mil. 

Ml?. 

Ml7. 

M  IM. 
Ml*. 
MIF,  . 

Ml  7. 
Mil. 
M19. 

M?  h • 

M  ?  I  • 

m:?  . 
M2I. 
•*?  M  • 
M2*. 
M2  F  . 

M  2  7  • 
M2  R. 

M  2  • 
m:*h. 

Mil  . 

M  72  • 
Mil. 
M1M  . 
Ml*-. 

M?r . 
Ml’’  . 

Ml  P. 
Ml°. 

M MU  . 
MM  l  . 
MM’  • 
MM7  . 
MMM  • 
MM*. 
MMH. 
MM  7. 
MMM  • 
MM  «>. 
4*0. 
Mf  l  . 
M*r. 

M  r  7  » 
M!  M  . 


-  TN  COMPRFSSTONI 


-  r*.  SHf  ARl 


rr  i ncho 1 1? i •  i?m •  ’ 

jf  f  K  FEI N01  »2S2*2Bf»2  6* 

5ftM  PRINT  287  »N0*  I  7Z«  MO  »  1 1  1 1=  1  *G  1 
287  FuR“ ATI tr  »C  c 1? 7KTFNST0WI 

10  TO  97  ,  e| 

-  c  c.  print  177  INfit  II  •  ■ 

; \ 8  r MM  7  •  K . sr  12 .« .8 . .  1 IHCOH-HF «I o* . 

r.O  TO  97  ,  .  _t 

r  PRINT  239*NO.<77INn»II«T-l*6l 
2  8®  FORMAT  1 TFtfcFl? . * . . *HR HF AR I 
00  TO  91 

12M  TF  IKFFINOI  I  92«SC#  EM 
3M  PRINT  87#NO 

C?  roftMAT  IEX.T7.13P  -  IN  TFNr  ION» 

00  TO  27 
pf  PkTNT  68 1  NO 
80  FORm*T  (PX»Il«l7H 
CO  TC  07 
02  PRINT  39*NO 
8®  rORH  AT  <  f  X  *  T  ?•  1 1  H 
T!  NNB-1 

TF  I  NTNN*2 -NLR  1M9.M9»M2 

L INN -NTNN 
0  0  MO  l-l  iM 

TF  IN0D1N0*  II-LINNIMl.Mli.  Ml- 
L INN-N001 M0«  T 1 
CGNT INUr 

NNO=«L  YNN-l  »  /NL*«  l 

IF  INNP-NBM* MM.M9.MM 
MIL  A -N3N-NN0. 

00  Ml  Jr  1 »  K  TL  A 
BAT  KSFAFP  Hi 
M  7  CUNT TNUF 

N2  =  NLR*? 

RFAOL0|l«l  <FINlf ISIN.K) fM=l .NflUl tN=Nl *N2I 
PftCKS PACT  m 

145  R  f  AO  ^P(  IT1  •  l,CFINI.ISlNtMI.M:l.NOKltK  =  l.NlPl 

T  1=  N^O  «  NO  tl  1 
|  2  r.NO  0<  NO*  3  1 

.«••!  RAO  I  T2  MRAll  Tl  I  1 

♦.  22,f*1RAOIT2I“RAO<I11  I 
12  T  T61  r21  -TTA  I  II  l» 

TP  1  12  I  29M  »29M  »29|S. 

?®M  y2:.**|P.2e7lO-TTAIIlll 
:nr.  y 1  =  T TM  C  Tl  I  ♦  I f 
F:PO«NOt  11 
PR=F0<  NO.I I 
CALL  TNT  PO 
no  mo  T  =  l»M 
r.o  to  7  orl.NNSO 

TF  <  NOOf NO  *  T I  -N0031 J1 IM7  »M€*MT 

to  f .  00  M r  K=l»? 

TF  I0ISP1  of*  I -IU».  I7N*M9  »M* 

7D  KP:  t  T-l  I  *2  ♦  K 
OU  71  L-l»R 


MI- 


MI 

MU 


MM 
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4*  •  . 
45  0. 
q«  7. 
« *;  3  • 
4»-9. 
Mo*:  • 
MCI  . 
4F?. 
4F7. 
qf  q. 
*er . 

qbn. 

qe 

43  3. 
«f r-. 
M7i:  • 
«l7I. 
47?. 
q  7*. 
q  7  q . 
47* . 
i7r. 
4  77. 
479. 
9  7  3. 
4  Q* I  • 
4  S  1  • 
4  b;. 

q  3  j  • 
494  . 

qer. 

q  3r.  • 

qp7. 

4  so. 

4  O':. 
4*>m. 

qr-i . 
4*3". 

1,07. 

494  • 
495. 


SI  IKP  .LI  -tl. 

SLi  L*KP  l  =  »i  • 

71  CONTINUF 
45  f  ONTTNUF 
47  CONIINUF 
43  CONTINUE 

r>o  ?s  K  =  i.q 

0  5  7  9  J  =  1  #2 
XN  =  (K-1»  *  2«  J 
KM= 1 NCOINO.K I -1  » • 2 ♦ J - 1 NN6  -1 l*NL6 
00  7*  L=l»4 

TF  I  NODI  NO  *L  l-qOUNOtKl  17  5*74  .74 

74  00  7r  MX  :  1*  2 

L  t  L -1  I  •  5  ♦  H  Jl 

L  K  =  I  NO !><  NO*  L  »  -NO C  I  *0 •  K  I  )  *  2  ♦  MX  -  J  ♦  l 
IF  I  L  M  I  7  e.  »7r  .7“ 

72  SfKPtLMl  rS(KM#LMI- *3*SL  tKN* IN  > 

75  CONTINUE 
7 J  CONTINUE 
7 *»  CGNTTNUF 

&ACK5F  ACC  I  I 

KM T  TF  (till  IFIN  l*IS<N,M|*M=|  «NBWl  *N=1»NL0» 

IF  I NNB-N  8N  *  32*  B  ?•  S2 

32  URITF  (111)  IF(M»ISIM.M|,Mri,*iawi  .NiMifN2l 
B ACKSP  ACF  H» 
no  175  K=l.MtLA 

REAP  (  lfi  I  I  F(N1  .  IS  «N*Ml  »M=t*NBW  I  «Nri  *NL&> 
l»e  CONTTNUF 
5?  NODINO.41  :-NOOINT»  41 
4  CO  FONTTMUF 
CALL  OISFL 
TF  I N'„  B  1 9  9  8  »  9  3  P  *430 
43li  CONTTNUF 
GO  TO  993 

"9?  PR  IN  T  2SF*NAMCHA 

2?F  FORMAT  <8M  FLFMPN-  T#T4  , 

1 4 1 1 H  WHICH  IS  TN  wrriMlTy  OF  LOAF*  HAS  FA  IL  FO/ 
lloH  RUN  TS  TFRMlNATFni 
956  STOF 
F  NO 


40  OF  COMPILATION 


NO  DIAGNOSTICS* 
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1 . 


*4  . 

-  • 
f  . 
7. 

P  . 

• 

IH. 
11  . 
1?  • 

I  *  . 

I  a. 

I*. 

IF-. 

17. 

1  S. 

II. 

?  o. 
21. 

2  2. 
21. 

2  ft  • 
2*. 
2r . 
?■». 

2  B  • 
21. 
?.n. 

ii. 
.*? . 
73. 
3ft  . 
3*.. 
If. 
37  . 

3  B  • 
3n. 
AO. 
<11. 
ft?  . 
ft  1. 
ftft  . 
1r  . 
ftf  . 
ft?  . 
«!  ?. . 
Ml. 
*  l  . 
51. 
3  2. 

/•  f 

..  .  • 

r  M  . 


SUBROUTINE  OtSPL 

COMPON/COLO  r/NBli  j  »  MM  »NU«BLK  »NF6  •  0  •  2til»>  •  A(  2Mt«  iOlil 
73- 1  ,ue-c 
MCNTG 
NF-  1 

NH-NN *NN 
MT=NF<1 
K.L  =NN  ♦  1 
NO-  1 

T P  (NUMpLK-l I Di» If  3* 01 
11  REMIND  11 
REMIND  111 

NS=M 

co  to  i*n 

r  reduce  equations  iv  blocks 

c  1.  SHIFT  nL  0  CK  OF  Fyu  AT  TOMS 

l(;f!  MS  -  NR*  1 

00  12f  N=l»NN 
NM=  NN* N 
B(Ni  :r(Nmi 
R  (  NM»=ll. 

00  l?r  M  -  1 1  MM 
A  (  N  »  M  I  -  A  f  NK  »M  I 
1  c*  A  « NM •  Ml  =(•  • 
c  2.  READ  not  block  into  CORF 

TF  (NUMRLK-NBI 1 fU»  Iff » l'U 
JK|>  READ  (1(11  (R  (N  I  .  (  A(N  *H  I  »s=  1  »MH  J  .NZML  »MH1 
TF  (KBIlf  Mt  lOltt  1  f  M 
C  7.  RFOUCF  BLOCK  OF  EQUATIONS 

If*  KT=NROK*?-(NUKFLK-il*NFQ 

NF=U 

ISM  00  3 eu  N= 1. NT 

C  CHECK  FOR  VF  R  i  SMALL  03  ZERO  ELEMENT  ON  DIAGONAL 

CHIU. 

0  0  loS  MF 1 t  MM 

IF  (FH-ABSI A (U • M  I  I 19M»lEf»lEE 
<)4  CH=ABS(A(NfM|  I 
Iff  CONTTNUF 

Xlk^AOSIAINtl  M/CH 
If  (X7X-7RI rft S5f 2C? 

C  TF  7rRQ  OR  SMALL  ON  DIAGONAL  TERMINATE  RUN 

7*  PRINT  2  T.Plif  XXXtNtNR 

•MllU  FORMAT  l™H  SMALL  DIAGONAL  -  A(N«1I/MA>  A(N.M»  ISt 
lE14,Otl2H  TN  FuUAT  TON.  Iftf  ?»!  Op  BLOCK.  TUI 
NF 

GO  TO  E(. 

C  NORMAL  ELIMINATION  PROCEDURE 

JOC  R ( N I  r  p ( N I  7  A ( N  »  1 1 
00  77S  L -  2  .MM 

I  P  (  ABS(  MK.LI  l-ZRI  3  75.  3  7ft  J7E 
T 7 f  C=  A  (  ML  1  7A  (  N»  1  » 

;=N*L-l 
J-  •'< 

DC  Tri>  K=L.PM 
J-J*  1 
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1UM. 
U.l . 


*£''  4<  t  »*;  i  =  a  ( r « j  i -c  *a in  »k  i 
ecn “Be  n-A(N»u*b'Ni 

A  *  N .  L  »  =  C 
I7C  CONTTNUF 
’oil  rOHTTHUC 

c  hsitf  block  op  FijuArroNS  on  tapf  n 

IP  fNUMRLK-NBl  f  M  •  4LiLi»  S«i 
^  ur<Trr  111  MR(NI.(AI(M.HI,M=1»HM|,N=1#NMI 
CO  TO  101) 

c  BACK -SUBSTITUTION 

“on  CO  AfU  b=1«MT 
N=MT«1  -M 
LI  =2 

97  00  U'5  K=L1«MM 
L  rN4K -  1 

,J  i  ?  0  (  N  I  =  R  (  N  I  -  A  (  N  ,  K  I  *  3 1  L  I 
BSO  CON'T  TNU F 

DO  HECl  N  =  1,KT 
NM=N*NN 
ACNK.MPJ  =p.(,\-| 

IP  l  NR-MUMRLK  l?l!tUfci],l|Cr» 

7 H  A  (NMtMFf =61  SMI 
<«r'l  R(NM|=6<N» 

N  8  -N  B  ~  1 

IF  (  Nfl  IBlltJ  .fijis,  71 
71  PACKSFACF  11 
MT=  NN 

FFAP  «  1 11  (  PCM  •  f  A  (  N  •  K  )  #  M  Z 1  *  *)  m  J  t  N  C  1  •  N  M 1 
RACKS  F'ACF  11 
CO  TO  ()ii|.i 
rtlll  K-  1 

NMI-NEu-1 
0  0  CMIJ  r=l#KROW 
NMIrNM  !♦  2 
NtfCr  NMT  ♦  1 

7F  C  M  K  F  -  ?  ,*MFm  |  £  ft  2  t  E02»  Fill 
SHI  K  =  K«l 

NMT=NCG* 1 
NMF=N«I«1 
fh?  roNTPriui 
Ffl 1 1  CONTINUE 
CO  RFTUON 
r  NO 

C’*\ . * . 

. . . 


£ of  compilation: 
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1  . 
t  . 
». 
V. 
r  . 
G. 
7. 
3. 
9. 
10. 
1  1  . 
1 

IT. 
1*. 
lr. 
ir.. 
17. 
1-1  • 
1  r. 
• 

?  i . 
2  2  . 
2*. 
?«l. 

<  T  • 
2  G  • 
2  7. 
23. 

<  p. 


GUP»0UTIN£  MAH. 

COMMON /JACX?/CM«*<  I  •  B  f  T#3I  »FL I  3  I  «  SL  I  8  *  8  I  *BP  I  2.81 
COKMDN/BCF  AP/FP.  f  .T.  TCP.  CC=f 
C  OMMGN/K  R  TNG / X 1  tX2.fl  •  *2  •  *  *  2  I  .  7(3  I  •  Zll*  • 

P.R=X  1  «X2*x  I  Ml 
njM=*»T/t 1 . -FR*FR I 
X  X  =X  <  M » 
fiz.il  N  I 

pi  i«  ? i =i./» ? 

(V  l  <4  UVV/X2 
Fit.  l»:i./Kfi 
HI  2 *2 l=XX /PR 
PIC. ? I ryy/RR 
n  2  «(i  irxx  .yy/Rp 
Mr. 7  *  - 1 ./RR/y  < 

212  »3l-XX/RR/y2 
PIT.  Tl  zn  c •  71 
M  T  .«  l:M  2  .9  I 
M2.  r  I  =-Ft  2.  H 
PI  1.G 1  =  1 . /X2 -XX/RS 
M?.  7I=-PI 2. ?l 
an .3 i=yy*Br* »f  i 

no  2F  T=l.S 

net i «r i=p jK*tn< i • t 2  «ii i 

PM?.  Ilrnjy.*  (nl2.CMFRvPll.nl 
961  3. 1  *11  .-PR  I  •«!  3.H 

2*-  CONTTNUP 
Rt  TURN 
FND 


£\I3  of  compilation: 
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1 . 


<4 . 
r . 

G. 

7. 

3. 

°. 

l'l. 

11. 

1*. 

17. 

1  **  . 
lr  . 
If. 

1  ■». 
13. 
lp. 
2'1. 

2  1  • 

2  7. 
27. 
2*. 

2  '  . 
2G. 

2  7. 
73. 

2  "  • 
711 . 

71  . 
72. 
;7. 
74  . 
?r  . 
7G  • 

’  7  • 
38  • 

7  «. 

qn. 

4  l . 
47  . 
4*. 
44  . 
4  r  . 
40  • 
4  7. 
43. 
4  n  • 
eli  . 
f  l  . 
e  2  • 
*  2  . 
*4  . 
*  f  . 


C0MM0N/JArKF/C(4.4 I.R<7f9l .TL181 * SLJ® * ® 1 ’0”‘ 1  *  8  * 

C0HMnN/Xft:NP/Xl«*2*n*  T*  *V 

OT^GNSTHM  HI3l.AAm.OU8l.ni3. 31.0)113. 31 

PAT  A  (HI  T».  T  =  i.  FI /.F*5f  E55E*. 8  8688888.. E5Ef - • 

OAT®  t*AII».I=l. 11/ -.775.0 •*•778/ 


PO  3  T-l.3 

Oil  1 1 -li  . 

FL  ( T  I  =li. 

0  0  4  K  =  1  .3 
ta  r>  (T.KI-l* » 
0  0  e  J=T.  .8 


CCt. Jl -n. 

*  SL( J.I  l=U. 

2  CONTINUF 
00  71  L  =  1  *  1 
K=L 

KCMI-AAIMI 
T?R  =)•  1*X'*X(M| 
00  7f  J  =  1  .3 


N  =  J 

y|  N  l  =  AA  IN  I 
CALI  HATF 

nUHl=H(Ll*H(J|.RR.X?.7? 
00  2  P  NROA=l»B 
0  0  20  NC0L=NR0H»3 


0UM2=li. 

no  in  rfHm-* 

OUM2  =  OUH2«B(KH.NRG4|.&MKK.nCOLI 

0(  NR0H.NC0L 1  =  0 INSO* .NC OL I  *  0UK2.&UM1 

CONT  TNUF 
CONTTNUG 
00  F  T  =  l»7 


J=  I  ♦  1 

00  e  K=W.3 
6  0 ( K  » T  l  =  0(  I  »*  I 
t  1  =-l 

0  0  2r  1  =  1  »8 
T  1  -  T 1  *  1 
T  2= ( T-l  1/2*1 
DO  24  J=  1 »  3 
00  27  K  =  1  »4 
K 1 =K  ♦  II*  4 

<7  0  X  (  T»JI  =  0MT»JMC(K»T2  I *0<  *1 »  J I 

24  CUNT TNUF 

TF  IT1-1I2C.  .21.21 


21 

?r: 


r  i  =~  i 

C ON T TNUF 

r  i  =- 1 

00  23  1=1.8 
1 1  =  r  i .  i 
T2  =  I T-l  1/5*1 
DO  27  o=  T . 6 
OU  2 G  K  =  1  * 4 

x :  =k  ♦  r  l  •  4 


u 
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f  r . 

2  T 

sl  ( j.  ii  rs:.  IJ.ri.nx{j.Kii*C<K.T2» 

r7. 

•5  *» 

.  « 

conttuuf 

f  r. 

T  r  t  Tl-l  Kf  #  22. 2  2 

7  ? 

T1  -  -! 

r  n. 

23 

COM  TNUF 

Ft. 

00  7  1=1,7 

f  2  . 

j:r»  1 

€  !. 

0  0  7  K  =  J,S 

f  *  . 

7 

5Lir.tOiSl.IK,  II 

FF. 

»F  TUPN 

f  r . 

CND 

or  compilation:  nc  diagnostics. 
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i\i  r*» 
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QUfiQOUTtNr  UMTFQMU.fi.fi 
nrR ANUN! R » 

C=A«  <8-6  I  *0 
RFTL’QN 
r  NO 
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